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Black strings, one class of higher dimensional analogues of black holes, were shown to be unstable 
, to long wavelength perturbations by Gregory and Laflamme in 1992, via a linear analysis. We 

revisit the problem through numerical solution of the full equations of motion, and focus on trying 
to determine the end-state of a perturbed, unstable black string. Our preliminary results show that 
5-H ' such a spacetime tends towards a solution resembling a sequence of spherical black holes connected 

by thin black strings, at least at intermediate times. However, our code fails then, primarily due to 
large gradients that develop in metric functions, as the coordinate system we use is not well adapted 
to the nature of the unfolding solution. We are thus unable to determine how close the solution 
we see is to the final end-state, though we do observe rich dynamical behavior of the system in the 
intermediate stages. 
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V/-^ ■ I. INTRODUCTION 

00 ■ 

' The stabihty of four-dimensional black holes is a well known, and fundamental result of relativity theory P|. 
. The picture in higher dimensional spacetimes was shown to be quite different by Gregory and Laflamme 0, Q , who 
' demonstrated the existence of unstable long wavelength modes of the black string in perturbation theory. This finding, 
, coupled to arguments based on entropy considerations, led to their conjecture that black strings might bifurcate into 
higher dimensional analogues of spherical black holes. Cosmic censorship would be violated were this the case, since 
^H, a singularity must be encountered by a bifurcating black hole horizon, essentially as a consequence of the principle of 
I ■ equivalence |4lj. 

^jpj' The existence of the Gregory-Laflamme instability has been assumed in many subsequent studies of higher dimen- 
• • . sional gravity theory, including the classical limit of string theory (see, for example, 0,0, Q,|Si])j some of which have 
. ^ ' ^^so assumed the validity of the bifurcation conjecture. However, a linearized analysis can say little, if anything, con- 
, cerning the nature of the full non-linear evolution of an unstable string, and the final end-state of such a configuration 
\^ ' remained to be established. 

5^ \ Recently, Horowitz and Maeda were able to prove, under some assumptions, that black strings cannot bifurcate in 
finite time 0- Furthermore, they conjectured that the system is likely to approach a new stationary solution which 
is not translationally invariant along the string direction. However, even if the assumptions involved in the proof 
are sufficiently generic, their analysis cannot identify the final end-state of evolution. Partial answers can be sought 
via perturbation analysis as done by Gubser |lQj |. By assuming the Horowitz-Maeda conjecture and linearizing the 
solution at the critical length to first order (and a partial extension to second order), Gubser argued that the transition 
to the final solution must be of second order type (i.e. discontinuous). Despite these developments, it seems clear 
that a convincing answer to the question at hand can only be obtained by solving the full equations governing the 
problem. A step in this direction would be to search for special solutions, such as stationary ones, and compare the 
physical content of the obtained configurations with the black string solutions. This has recently been carried out by 
Wiseman 11], who numerically solves the equations resulting from a static ansatz. Interestingly, he finds non- uniform 
solutions with mass larger than that of the black string for a given compactification radius. Thus, he concludes that 
the solutions he finds cannot be the end-states conjectured by Horowitz and Maeda (also see related work by Kol |1^|). 

Additional work by Unruh and Wald [13| studies the dynamics of a uniform cylindrical configuration of matter in 
Newtonian gravity. They observe that a perturbation of the density gives raise to a Jeans instability responsible for 
the collapse of the system along the cylinder's length. They then argue that if the main features of this model are 
robust in the passage to the general relativistic system, one possible end-state for the perturbed black string would be 
collapse in the string direction, resulting in singularity formation. Note that this collapse need not lead to violations 
of cosmic censorship, as the final singularity could still be hidden by an event horizon [l^ . 
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Clearly there are several distinct viable possibilities for the final end-state of a perturbed black string, with remark- 
ably different consequences associated with the range of options. Current conjectures range from "nothing interesting 
happens" , to violations of cosmic censorship, to the arguably more extreme case of a complete collapse of the space- 
time. In order to completely settle the issue, the full dynamics of the perturbed black string needs to be addressed. 
At least in principle, this will allow us to identify which of the above possibilities (if any!) is actually realized. In 
this paper we report on preliminary work in this direction — a program to simulate the dynamics of the black string 
through numerical solution of the Einstein equations. At this stage of the project, we cannot yet provide an answer 
to the question of the end-state; however, we have tantalizing results that show the spacetime going through a config- 
uration resembling 5-dimensional spherical black holes connected by thin black strings that expand along the string 
dimension. Our simulations eventually crash while the spacetime is still fairly dynamical, and so we cannot determine 
whether what we see is near the end-state, or merely an intermediate configuration in a more complicated evolution. 
Underlying the current failure of our simulations is the fact that the coordinate system employed is not well adapted 
to the solution that unfolds at late times, wherein fatally steep gradients develop in metric functions ^. 

The outline of the remainder of the paper is as follows. In Sec. m we begin by describing the equations of motion, 
our coordinate choices, generation of initial data, as well as our numerical solution scheme. Additionally, we also 
briefly mention the tools we employ to monitor the solution, deferring details to the Appendices. In Sec. IIIII we 
discuss the results obtained with this code, and conclude in Sec. IIVI bv mentioning directions for future work that 
may allow us to more definitively answer questions regarding the end-state of the Gregory-Laflammc instability. 



II. EQUATIONS, BLACK STRINGS AND NUMERICS 

We wish to solve the vacuum Einstein equations in higher dimensional settings. For simplicity, and without loss 
of generality in studying the Gregory-Laflamme instability, we only consider the 5-dimensional case, and restrict 
attention to spherical symmetry within the 4-dimensional subspace tangent to the "extra" dimension. We also use 
the natural generalization of the ADM decomposition to derive the system of equations that we then solve numerically. 
Choosing units in which G = c = 1, and adopting MTW |l5j| conventions, our starting point is thus a metric element 
given by 

ds^ = {-a^ + -iABl3^fi^)dt' + 2-fABf3^dx^dt + jAsdx'^dx^ + -fndn^ (1) 

where — (r, z), and dV,^ is the 2-spherical line element with coordinates chosen orthogonal to the t — constant 
congruences (hence there is no shift corresponding to angular directions). All metric components defined via 
depend upon {t,r,z): t is a time- like coordinate, r is a radial coordinate, and z is the coordinate along the length of 
the string. To further expedite the numerical implementation, we make z a periodic coordinate by identifying z = 
and z — L. Then, following the results of Gregory and Laflamme, we expect black strings to be unstable only for L 
greater than some critical length L^- 

The vacuum Einstein equations, written in ADM form 0,0], are 1) the Hamiltonian constraint 



2) the momentum constraints 

3) the evolution equations for the ^ab 



H ='^'^^ R + K^- KabK"'' = 0, (2) 
Ma ^ Kj' I, ~ if |„ = 0, (3) 



b 



dt 



-2aKab + Pa\b + I3b\a, (4) 



that follow from the definition of the extrinsic curvature Kab associated with t — const, slices, and 4) the evolution 
equations for the extrinsic curvature 



+/3Vifeb + P''\bKc.a + l3'K,b\c + aF^F^j^dH. (5) 



^ This divergence of metric gradients does not happen earher as resolution is increased, and is not accompanied by divergence of curvature 
invariants such as the Kretschmann scalar. This suggests that the code is evolving to a coordinate, rather than a geometric, singularity. 
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In the above, a,b, . . . are four-dimensional (spatial) indices, '^^^ Rab and ^^^i? are, respectively, the Ricci tensor and 
Ricci scalar intrinsic to the four-dimensional spatial hypersurfaces, a is the lapse function, (3'^ is the shift vector, the 
vertical bar | denotes covariant differentiation in the spatial hypersurfaces (compatible with 7ab), and = —2S!^S^. 
We note that the term proportional to iJ in ||SJ) has been added as a result of stability considerations; see for instance, 
the discussions in 17, 18, 19.] . To simplify the final set of equations solved numerically, as well as to regularize certain 
terms that otherwise diverge at spatial infinity (see the discussion of our coordinate system in the next sub-section) , 
we define the following variables 

9rr = T^r 9rz — ^rz Qzz — ^zz 

gee = m/r"^ 94,4, = m/ir^ sin^ 9) 

krr = ^ Krrj k^z = ^rz kzz = -^zz 

kee = Kgg/a k^,^, = Kgg/ {a sin^ 6) (6) 

and use them as the fundamental dynamical quantities in our numerical code. As discussed in the following sub- 
sections, to complete the prescription of the evolution problem we need to choose a suitable lapse and shift, specify 
initial and boundary conditions, and then implement these choices and specifications consistently. 

A. Boundary and Coordinate Conditions 

A particular concern here is that "standard" outer boundary conditions often imposed during numerical 

evolution of Einstein's equations, might not be well suited for studying the string instability. In particular, we must 
be able to evolve for very long times while absolutely minimizing spurious influences from the outer boundary of the 
computational domain. In addition, in the present case we cannot assume, a priori^ that any given initial configuration 
will settle down to some stationary solution; thus, boundary conditions predicated on such assumptions (such as a 1/r 
fall-off condition in a metric component), when imposed at a finite proper distance from the black string, could very 
well adversely affect the numerical results ^. To ensure minimal boundary influence we therefore extend the domain 
of integration to i° by radially compactifying the spacelike hypersurfaces via the introduction of a new coordinate, x, 
deflned by 

r 



As might be expected, this transformation causes computational problems of its own — most notably decreased spatial 
resolution at large distances — but, as discussed in Sec. Ill CI we can deal with these difficulties using numerical 
dissipation. Having introduced the new compactifying coordinate, we can directly impose boundary conditions derived 
from the demand of asymptotic flatness at spatial inflnity, which lies at a; = 1. 

We employ singularity excision techniques [2iJ| to allow us to evolve the entire perturbed black string spacetime 
exterior to the apparent horizon (plus a certain "buffer zone" that lies within the horizon). Hence, we do not 
need to impose inner boundary conditions as long as the t = const, hypersurfaces penetrate the horizon, and that all 
characteristics of the evolution equations are in-going on the boundary. Ensuring that this is the case involves choosing 
"good" coordinate conditions (choice of lapse and shift), which, for generic string evolutions, remains an open problem. 
As a preliminary step, we have based our coordinate choices on those that yield the ingoing Eddington-Finkelstein 
form of the unperturbed black string metric 

dsls = -(1 - 2M/r)dt^ + Ul/rdrdt + (1 + 2M/r)dr'^ + dz^ + r^df]^ ^ (g) 



^ In fact, in an earlier version of the code that did not use a radially compactificd coordinate system, we did encounter such problems, 
in that some artificial stationary non-homogeneous solution was apparently entirely "sourced" via an outer boundary located at a finite 
distance from the string. 
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Comparison with the general 5-dimensional ADM form provides the identifications 

aBS = (1 + 2Af/r)(-i/2) , (9) 
/3^S = i2M/ir + 2M))6^. (10) 

For reference we also list the two non-trivial components of the extrinsic curvature of a.t — const, slice defined by ((HJ: 



r + 2M ' 



In generalizing jnj and to the dynamical case, we have chosen a = ass and — /9gg = 0. In a preliminary 
version of our code, we also required that /?'' = /Jgg. This, however, caused a coordinate pathology to develop at late 
times during the evolution of unstable strings — specifically, some regions of the horizons approached a zero coordinate- 
radius, while maintaining finite proper radius. In our current efforts, we choose (3^ such that ggg remains constant 
during evolution ,2^ , by requiring that 

"tee r 



This shift condition performs reasonably well, as will be seen in Sec. IIIII However, our current simulations still suffer 
from "grid-stretching" problems in the z direction at late times, suggesting that a more dynamical gauge condition 
for (and possibly for a and (3^) could be useful. This issue is discussed in more detail in Sec. II VI 



B. Initial Data 



As anticipated, we observe that even numerical truncation errors, if non-uniform in the z direction, are enough to 
trigger the Gregory-Laflamme instability in our simulations. However, to reduce the computational effort required to 
reach the "interesting" (i.e. non-perturbative) stages of evolution, we adopt initial configurations whose departure from 
the black string solution can be arbitrarily tuned. In order to find such data, we must solve the Hamiltonian constraint, 
and the r and z components of the momentum constraint (the other components of the momentum constraint are 
trivially satisfied because our coordinate system is adapted to spherical symmetry). The deviation — not necessarily 
small — from the black string solution, is introduced via ggg, and takes the following form: 

ffee(0,r,z) = l + Asin(^z^^e-("°)'/^'. (13) 

Here, A is used to control the overall strength of the "perturbation" , while q is an integer that controls the spatial 
frequency in the z direction. For the results presented below, A = 0.1, q — 1, tq — 2.5 and dr — 0.5, and we perturb 
about a unit mass (Af — 1) black string solution. As described in more detail in Appendix 1X1 g^r, krr, kgg are then 
calculated by solving the constraint equations, with the remainder of the metric and extrinsic curvature variables set 
to the values they would take for an unperturbed black string (see © and (fTT|l ). 



C. Numerical Evolution 



To numerically evolve the initial data sets described above, we discretize the evolution equations Q-lO using 
second-order accurate finite difference techniques that include Crank-Nicholson treatment of the temporal and spatial 
derivatives. We use a uniform distribution of grid points in z and x (recall that x = r/{l + r)). The resulting implicit 
system of algebraic equations is solved iteratively. We initially implemented a serial version of the algorithm, and later 
coded a parallel version using the CACTUS Computational Toolkit 23] , wherein the equations of motion, monitoring 
tools and I/O were handled by our own routines, suitably interfaced to CACTUS. 

Black hole excision is handled as follows. We periodically find the apparent horizon, as discussed in the following 
sub-section and Appendix^ We then define the surface along which we excise to be a certain number of "buffer" 
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points inside the apparent horizon (typically 10 — 30 buffer points are used) ^. During each Crank-Nicholson iteration, 
all the evolution difference equations are applied up to the excision surface, and any function values referenced by 
finite difference stencils interior to this surface are defined via fourth order extrapolation. When the apparent horizon 
location, and hence excision surface changes during evolution, function values at all repopulated points (i.e. those 
that moved from inside to outside the excision surface during the time step) are computed via the same fourth order 
extrapolation routine. The one exception to this procedure is for the grid values of ggg, which we specify a priori on 
the entire computational domain, and that remain fixed due to our gauge choice (|12|l . Moreover, we have found it 
useful to choose a functional form for ggg that tends to zero at some positive value of r (though inside the original 
apparent horizon location and outside of the limits of integration of the initial data) . This causes the "pinching-off" of 
the unstable black string to be less severe in coordinate space, i.e. we approach zero areal radius at a finite coordinate 
r. In turn, this slightly reduces the virulence of the coordinate problems we observe at late times, and also provides 
better load-balancing of the parallel code, given the method CACTUS uses to distribute grids among processors. 

For the evolution, we choose a time step dt = AcFLniin(dr, dz), where the constant Acfl must be set less than 
1/V2 in order to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition that results from our iterative solution 
of the Crank- Nicholson scheme (typically we use Acfl = 0.25). Note that this restriction on Acfl is based upon 
flat-space light speeds within our coordinate system (0 with M — 0), which, for the solutions presented here, are 
always greater than or equal to the actual coordinate light speeds. The function min(dr, dz) is calculated by only 
considering mesh spacings within the non-excised portion of the coordinate domain. Thus, as the excision surface 
moves, dt changes with time since our grid is uniform in x, and hence non-uniform in r. 

Crucially, we add Kreiss-Oliger-style [23| numerical dissipation to the evolution equations to control unphysical 
high-frequency solution components ("noise") that would otherwise arise during the simulations. This is particularly 
helpful at the excision surface, and near i'^, where the radial compactification of points causes all outgoing wave-like 
components of variables to eventually become poorly resolved. Smoothing of the high frequency components via the 
Kreiss-Oliger dissipation — which only targets wavelengths of size on the order of the mesh spacing — prevents them 
from inducing numerical instabilities near the outer boundary. 



1. Monitoring the Evolution 

To elucidate the nature of our computed spacetimes, we monitor the following quantities: 1) the location of the 
apparent horizon (which is also used for excision as explained earlier), 2) the trajectories of null geodesies that, to 
a certain extent, should trace the event horizon, and 3) the Kretschmann invariant / (the square of the Riemann 
tensor) : 

/ = Ro.p^sR'"'"''. (14) 

If cosmic censorship holds — and results from our current simulations provide no evidence to the contrary — then any 
apparent horizon found will always be inside an event horizon. As is well known, although the apparent horizon can 
often be used as a reasonable approximation to the event horizon, the two do not, in general, coincide Clearly, the 
event horizon is the quantity of interest in studying the Gregory-Laflamme instability, and therefore we would like to 
locate it, or at least a reasonably good approximation to it, in our simulation results. Such an approximation can be 
obtained by looking for the boundary of the causal past of some r = const, surface that is sufficiently far outside the 
apparent horizon that it is certain not to be inside the event horizon, yet close enough to the apparent horizon that 
its causal past, tracing backwards from the end-time of the simulation, probes the region of interest of the spacetime. 
We use a method to find the approximate event horizon discussed by Libson et al. The approach is based on 

radial outgoing null geodesies; as explained in [25j |. the stable direction for the integration of null rays that emanate 
from the vicinity of an event horizon is backwards in time t. Thus, once we have the complete data from the entire 
evolution, we start with data from the latest time step available, and trace the null rays backwards in time. 

Monitoring curvature scalars is useful in obtaining coordinate independent information about a numerical solution 
of the Einstein equations. In particular, /, as defined by (|14|l evaluates to Ibs — '^^M'^ /■^gg^ for the unperturbed black 
string solution, and as jgg is an invariant in spherical symmetry (i.e. it is proportional to the area of an r = const. 
2-sphere), we can compare Ibs to the values of / computed from a numerical solution to get some indication of how 



^ In several tests, we also adopted an excision region given by the global minimum r value of the apparent horizon and compared the 
results with those obtained when the excision region was defined by the apparent horizon. The agreement obtained gives extra indication 
that the excision implementation is consistent. 

* Indeed, depending on the slicing an apparent horizon need not exist at all l26l . 
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close the computed solution is to a black string spacetime. Furthermore, we can examine / to see whether curvature 
singularities (other than the central r — singularity) may be forming prior to the demise of the simulation that 
invariably occurs when sufficiently steep metric gradients develop. We note, however, that if / does not diverge, it 
does not necessarily follow that the geometry is remaining non-singular; we would need to examine a larger set of 
curvature scalars to be certain that the solutions are remaining free of physical singularities. 

AppendixElcontains an explanation of the method we used to find apparent horizons, while details of the integration 
techniques aimed at approximately locating event horizons can be found in Appendix lO 

III. RESULTS 

In this section we present results from our preliminary study of the black string instability. After briefly showing 
that we recover some of the key Gregory-Laflamme results in the next sub-section, we present a detailed analysis of 
a typical unstable case in Sec. IIII Bl In the following, we will use the value of Lc ~ 14. 3M (with M — 1) found by 
Gubser [To| . which is more accurate than the value we can estimate from the zero-crossing of the (positive mode) 
interpolating curve presented in 0. 

A. Recovery of Gregory-Laflamme Results 

We ran a variety of simulations of black strings that were perturbed according to the prescription discussed in 
section III Bl We concentrated on cases with L ranging from Q.GLc to l.SLc and defined gee via equation (|13|l . In 
general, we observed the expected instability for L > L^, though for the maximum resolution at which we performed 
this survey (800 grid points in r and 200 points in z), we could only confirm Lc to within about 2% of the expected 
value. In this regard we note that as L approaches Lc from above, the growth rate of the instability goes to 0, 
requiring longer time integrations to identify the instability, which, in turn, demands ever increasing resolution to 
counter the effects of accumulating numerical errors. Furthermore, the initial configurations we have adopted contain 
energy in the form of gravitational waves, and some of this energy falls into the string early on during the evolution. 
The increase in the mass of the string (based upon the increase in area of the apparent horizon) is typically around 
0.3 — 0.5%, and we would have to take this into account were we to attempt to determine Lc from our simulations to 
a higher degree of accuracy. 

As a demonstration of the ability of our code to "bracket" the instability, and following the notation of ^3 : Fig. 
shows a plot of, A, defined by 

for L — 0.975Lc and L — l.OSLc- In the above, i?niax and -Rmin are the maximum and minimum areal radii, 
respectively, of the apparent horizon at some t ~ const, slice of the spacetime. In particular, we have A = for the 
static black string spacetime. 

B. Beyond the Linear Regime 

We now present more detailed results from the simulation of an unstable black string evolution. Specifically, we 
take L = lALc, since it is expected that this particular range for the z coordinate will yield something close to the 
fastest growth rate for the shortest wavelength instability Because we are now probing uncharted territory with 
our computations, we rely on convergence tests to provide an intrinsic measure of the level of error in our calculations. 
To that end, we ran the simulation at several resolutions {rir y- n^): 200 x 50 (grid spacing /i), 400 x 100 {h/2), 
800 X 200 (h/A), and 1600 x 400 {h/8). Due to our use of a compactified radial coordinate, the lowest resolution 
calculation cannot adequately resolve the late-time behavior of the solution. However, for the "medium resolution" 
computation with mesh spacing /i/4, we are apparently within the convergent regime — see Fig. |21 below for plots of 
the maximum and minimum areal radii of the apparent horizon as a function of time, as well as the quantity A defined 
by H15() . and Fig. Ofor plots of the norm of the Hamiltonian constraint as a function of resolution. Therefore, unless 
otherwise noted, all the results shown below are taken from the h/A simulation. 

Fig. 01 shows embedding diagrams of the apparent horizon at several times during evolution of the string, and Fig. 
|5l shows the proper length of one period of the apparent horizon (suppressing the angular coordinates) versus time. 
(Our embedding uses the vertical axis to represent the areal radius of the apparent horizon — the horizontal axis is 
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FIG. 1: The maximum (ilmax) and minimum (ilmin) areal radii, and the corresponding function A of the apparent horizon as 
a function of time, from the evolution of perturbed black strings with L = l.OSLc and L = 0.975Lc. The initial fluctuation 
in the plots correspond to the effect of the initial gravitational wave perturbation, most of which either falls into the string, 
or escapes to infinity. This close to the threshold Lc, the growth/decay of the remnant perturbation is quite slow, and so we 
cannot feasibly (at the resolution of the these simulations — 800 x 200 points in r x z) follow the evolution for much further 
than shown while maintaining reasonable accuracy (though we see no signs of numerical instabilities in the stable case, and 
such simulations have boon followed to 10, DOOM). However, the main purpose of this figure is to demonstrate the qualitative 
recovery of the expected threshold behavior for the onset of the instability at L = Lc- 



then uniquely determined by requiring that the length of the curve be equal to the proper length of the horizon) . The 
simulation crashes shortly after the last time frame shown, apparently due to the coordinate pathologies that have 
been discussed previously. The embedding diagrams suggest that, at least in the vicinity of the apparent horizon, the 
solution is tending towards a spacetime that can be described as a sequence of spherical black holes connected by thin 
black strings. Additional, quantitative, evidence for this conjecture can be obtained through a computation of the 
curvature invariant, /, on the apparent horizon. For an exact black string solution, this quantity, which we denote 



-^BS' 1^ 



/^s = (16) 
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FIG. 2: The maximum (7?max) and minimum (Rmin) areal radii, and the corresponding function A = (-Rmax/-Rmin — l)/2, of the 
apparent horizon, as a function of time, from the evolution of a perturbed black string with L = lALc- h labels grid spacing; 
hence smaller h corresponds to higher resolution. This plot, combined with the results shown in Fig. |H] suggest that the code 
is in the convergent regime — in particular at later times — for the h/2 and higher resolution simulations. 



while for the 5-dinicnsional spherical black hole, the equivalent quantity, I^n, is 



/^H = (17) 



AH 



where, in both of the above expressions, -Rah is the areal radius of the horizon. In Fig. we plot the normalized 
quantity 



evaluated on the apparent horizon of our numerical solution of the unstable spacetime — I^^ is 1 for a black string, 
and 6 for a black hole. The figure shows that, as judged by /°^, the part of the apparent horizon that is forming 
a long neck always resembles a black string — the part that is forming a bulge, however, has a value of I^'^ tending 
towards that corresponding to a black hole. At Rmax, I^'^ has only reached ~ 5 by the time the simulation ends; 
however, the behavior of i?max seen in Fig. |21suggests that the growth in the normalized curvature invariant, though 
slowing down, should continue. Fig. [Sjalso demonstrates the grid-stretching problems that we surmise are causing the 
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FIG. 3; The logarithm of the ^2-norm of the Hamiltonian constraint as a function of time, evaluated on the portion of the 
computational domain lying exterior to the apparent horizon, and from simulations at several resolutions of a perturbed black 
string with L = lALc- As with Fig. ^ this plot provides evidence that convergence is quite good for the h/2 and higher 
resolution simulations (at least until very close to when the supposed coordinate singularity forms, near t = 165). 

code to eventually crash — in that plot we use the coordinate z as the horizontal axis, and observe that the relatively 
small region where /"^ « 1 corresponds to the long neck in Fig. ^ In particular, in the vicinity of the "neck", gzz 
becomes quite large, as do its derivatives. 

Finally, in Fig. [7|we show plots of the approximate event horizon (as described in Sec. Ill C together with the 
apparent horizon for the simulation. The results shown in the plot suggest that our computed apparent horizon is an 
excellent approximation to the event horizon, at least at early times (not much can be said regarding the late time 
behavior of the event horizon, as the spacetime has not settled down to a stationary state when the simulation ends). 

IV. CONCLUSIONS 

We have performed a preliminary numerical study of the instability of the 5-dimensional black string. Coordinate 
pathologies prevent us from definitively identifying the final end-state(s) of an unstable black string. This claim is 
supported by the fact that the code crashes at very nearly the same time at varying resolution, and that curvature 
invariants remain well behaved throughout the evolution. The former suggests that a numerical instability is not 
responsible for the crash, while the latter indicates that a physical singularity is probably also not to blame. Despite 
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the premature termination of the simulation, we find evidence that the spacetime evolves towards a configuration that 
looks like a sequence of black holes connected by thin black strings, and characterized by an expansion of the string 
direction. Since the spacetime is still fairly dynamical at the time our simulations end, we cannot deduce how close 
this state is to a final configuration. Nevertheless, the dynamical behavior observed is sufficiently robust for some 
comments to be made. For instance, the results are not inconsistent with Gregory and Laflamme's conjecture that 
the solution bifurcates into a sequence of black holes — indeed, we can suggest at least two mechanisms by which this 
could occur: 

1. Via a thinning neck that eventually vanishes, if the trend seen in the simulation continues. Note that this would 
require a) that the proper length of the string continue to grow, in order not to violate area theorems, and 
b) that the thin string be non-uniform, for otherwise it would be subject to a further Gregory-Laflamme-like 
instability. 

2. Via a sequence of Gregory-Laflamme instabilities, if the thin neck stays "close" to a uniform black string, since 
the neck's length is beyond the critical one for a string with an effective mass computed from the radius of the 
apparent horizon. In this case, one could envision a "cascade" of instabilities leading to the bifurcation. 

We note that either scenario would not necessarily be inconsistent with Horowitz and Maeda's results, should the 
vanishing of the neck take infinite affine time as measured by local null generators of the horizon. At the same time, 
a continuation of the observed trend would argue against achieving a stationary solution with a mild dependence on 
the string dimension (i.e. a small value of A), as found in perturbative calculations For then, the rather extreme 
thinning/bulging that we see must be transient behavior that is "further" from the end-state than the perturbed black 
string was. 

A more complete exposition of the nature of unstable black string evolution would appear to require coordinate 
conditions able to adapt to solution features that develop at late times — that is, in a manner that does not introduce 
severe metric gradients that are not correlated with large gradients in physical quantities. For example, it may help 
to replace the fixed-lapse slicing with maximal slicing, which enforces that the divergence of the local, spatial volume 
element be zero. Another, perhaps even more crucial option, would be to introduce a z component to the shift vector 
that keeps gzz close to (or exactly) unity throughout the evolution. These options are currently under investigation. 

Additionally, it would be interesting to explore a wider range of initial conditions describing "perturbed" black 
strings than that considered here. For example, the imposed z-periodicity implies that the equivalent, uncompactified 
space time consists of identical spherical-black-hole/black-string segments at late (intermediate) times. It would be 
instructive to see what happens should we break this symmetry, by making L ^ Lc, and then introducing some higher- 
wavelength perturbation similar to q > 1 in H13|) . but with more asymmetry in the initial data (note that this would 
be more computationally demanding). Finally, it would be very interesting to study the evolution of the solutions 
recently found by Wiseman mentioned in the introduction. These configurations actually correspond to 

stationary solutions, and their perturbative stability, or otherwise, is currently not known. Since Wiseman shows 
that his solutions cannot be the end-states conjectured by Horowitz and Maeda, it is important to understand their 
behavior, since if they are stable they may well represent physically meaningful states, while if unstable, they may 
be difficult to attain via dynamical evolution. An interesting observation from our simulations, to the length so far 
achieved, is that they do not display a conical "waist" like those presented in and further analyzed in j27| . 

V. ACKNOWLEDGMENTS 

We gratefully acknowledge support from the following agencies, institutes and grants: NSERC, NSF PHY-0099568, 
The Canadian Institute for Advanced Research, The Canadian Institute for Theoretical Astrophysics, The Pacific 
Institute for Mathematical Sciences, The Government of the Basque Country, The Izaak Walton Killam Fund and 
Caltech's Richard Chase Tolman Fund. Computations were performed on (i) the vn .physics .ubc . ca cluster which 
was funded by the Canadian Foundation for Innovation (CFI) and the BC Knowledge Development Fund; (ii) LosLobos 
at Albuquerque High Performance Computing Center (iii) The high-performance computing facilities within LSU's 
Center for Applied Information Technology and Learning, which is funded through Louisiana legislative appropriations, 
and (iv) The MACI cluster at the University of Calgary, which is funded by the Universities of Alberta, Calgary, 
Lethbridge and Manitoba, and by C3.ca, the Netera Alliance, CANARIE, the Alberta Science and Research Authority, 
and the CFI. We would like to thank G. Horowitz, W.G. Unruh, R. Wald, R. Myers, T. Wiseman and B. Kol for 
stimulating discussions. 



11 



APPENDIX A: INITIAL DATA SOLVER 

We solve the set of coupled constraint equations via an iterative procedure, where at each sub-step of the 

iteration we solve a single equation for one of g„, kgg or k^r, assuming that the values of the other variables are 
known. We iterate this process until the residuals of all the equations are simultaneously below a certain tolerance — a 
typical value is 10~^. The overall iteration is initialized using values corresponding to an unperturbed black string 
solution. We now provide a few more details concerning the solution of each of the constraint equations. 

The equations are discretized on a uniform grid of points {xi, Zj) with i = 1, N^, and j = 1, Nz (recall that, 
from {T)), the radial coordinate, r, is related to a; by r = x/{l — x)). We first consider the Hamiltonian constraint Q 
which, in the coordinate system we have adopted, has the following form: 



dgrr , ^ d^9rr , Ogrr , f 99'. 



F,l^ + F2grr^^ + Fsgrr^ + Fi j + F5 ^^effrr = . (Al) 

Here, the Fm, m = 1, ... ,6, are functions that generally depend on all the metric coefficient and their derivatives 
except grr (and its derivatives). We discretize this equation to second order in the mesh spacing using a difference 
approximation centered at the points {xij^i/2i Zj). Because of the discretization used and the form of ljAl|l . we can 
solve the resulting set of equations "line- by-line" in x, starting at the inner boundary, i = 1, which is chosen well 
within the horizon of the string (typically at r = A/), and where the boundary values, [5rr]ijj j = 1; • • -^z, are 
those corresponding to an unperturbed black string. As we integrate outwards in x, the determination of the z-th 
line of unknowns, [grr\i,ji J = 1, ■ • -N^, involves the solution of a TV^-dimensional, non-linear, cyclic (because of the 
2;-periodicity) , tridiagonal system that we solve using Newton's method and a cyclic tridiagonal linear solver p^ . 
We now direct attention to the r momentum constraint, which, viewed as an equation for kgg^ has the form: 

Gi^ + G2fce9 + G3 = 0, (A2) 

where the G™, m = 1,2,3 do not depend on kgg or its derivatives. We note that there is complete decoupling in 
the ^-direction in this case; in effect, we have to solve an ODE along each z — const, line. We again discretize using 
second-order finite difference techniques, fix the boundary values [kge\N^,j, j = 1,----/Vz, at a; = 1 (io) using the 
unperturbed black string solution, then solve for the remaining unknowns, marching inwards in x. 
Finally, the z component of the momentum constraint, which fixes fc^r, has the form: 

dk 

Hi—!^ + H2k,.r + H3 = 0; (A3) 
oz 

where the Hm, tti — 1,2,3 are independent of krr and its derivatives. This equation is solved analogously to the 
r momentum constraint, but now using a discretization that is centered at points (xi, Zj 4.1/2 )• "Boundary values", 
[fcrr]i,ij 1 = 1,.. .Nx, are specified along the line z — Zmin, again using corresponding values from the black string 
solution, and the integration proceeds for = 2, 3, . . . Nz- 

APPENDIX B: FINDING APPARENT HORIZONS 

We use a flow, or level-set method to search for apparent horizons within t = const, spatial slices of the spacetime. 
We restrict our search to simply connected apparent horizons that are periodic in z. Such an apparent horizon can 
be described by a curve in the (r, z) plane, which we define to be the level surface = of the function 

F{r,z) = r - R{z). (Bl) 

In other words, the apparent horizon will be given by the curve r — R{z). The apparent horizon is the outermost, 
marginally trapped surface; hence, we want to find an equation for R{z) such that the outward null expansion, normal 
to the corresponding surface F = 0, is zero. To this end, we first construct the unit spatial vector s'', normal to 
F = const.: 



ab p 



(B2) 



Then, using s° and the t = const, hypersurface normal vector n", we can construct future-pointing outgoing(-l-) and 
ingoing(— ) null vectors 

r± = n'^ ± s'^. (B3) 
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The normalization of the null vectors is (arbitrarily) t'^l^a = —2. The outward null expansion 0+ is then the 
divergence of projected onto an F = const, surface: 

9+ = {g"'' - s-^s') Vd+a. (B4) 

Substituting expressions (|Bip and ljB2|) into equation ljB4p . with 6*+ = 0, provides us with an ordinary differential 
equation for R{z). 

Initially we solved equation ljB4p via a "shooting" method — given a guess for i?(0), and assuming that -R'(O) = 
(where a prime denotes differentiation with respect to z), we integrate the equation to z = z^, and repeat the process 
until we find a solution where R{zl) = R{z) and R'{zl) = 0. An efficient sequence of guesses can be generated using 
a bisection search, as the qualitative behavior of the solution is different depending upon whether the initial guess for 
i?(0) is inside or outside the apparent horizon. 

Since the shooting method is difficult to extend to a parallel implementation in an efficient way, we opted to use the 
following point-wise relaxation method (or flow method) to determine i?(z). We supply an initial guess, Ro{z), for the 
entire function R{z), and then iterate the following equation until the norm of the expansion 0+(z) of i^(0) is below 
some desired threshold (in our runs we have typically set it to lO^^h, where h is the basic scale of discretization): 

ARjz) = Rn+i{z) - Rn{z) = -0+{z)At. (B5) 

Here, Rn{z) is the solution after the n^^ iteration. Equation ljB5|) can be viewed as an explicit discretization of a 
parabolic evolution equation for R{z,t), where r is "time" and At is the time-step for the evolution (the parabolic 
nature of the equation is evident when 9^ is expanded in terms of R{z) via HB4|I — for brevity we do not give the 
explicit form of the equation here). Thus, for stability At must be chosen to be less than Az^. 

Given a "reasonable" initial guess Ro{z), one can see how iteration of equation (|ij5|l will cause Rn{z) to "flow" to 
the apparent horizon: if Rn{z) is outside of the apparent horizon, then typically the expansion 0+{z) will be positive 
there, causing i?„+i(z) to decrease towards the apparent horizon, and vice- versa if i?n(z) is inside the apparent 
horizon. We use R{z) = 2 as the initial guess at t = 0; after t = we search for the apparent horizon every N time 
steps (where A'' is typically in the range 10 — 30), and use the shape found at the previous search as the initial guess 
for the next search. Usually, on the order of tens to thousands of iterations of IjBSp are required to solve for R{z) to 
within a level of accuracy such that the approximate solution is roughly within a mesh spacing of the exact solution 
(as estimated in a few specific calculations by solving (jB5|l close to machine precision). A single iteration of ljB5|) 
can be computed very rapidly relative to the time taken to compute a metric-evolution step; however in a parallel 
environment, if thousands of iterations are needed on a regular basis (which is so at late times during the evolution 
of an unstable black string), the apparent horizon finder becomes a slight speed bottleneck in the code, due to the 
time it takes to communicate the results of each iteration amongst the processors involved. 



APPENDIX C: FINDING (APPROXIMATE) EVENT HORIZONS 

Here we describe one method we use, following |25l |. to locate approximations to event horizons. This method 
involves locating the boundary of the causal past of some r = const, surface of the spacetime by following radial null 
geodesies. 

We write the geodesic equation in Lagrangian form: 

>C = 9tt {t' f + 2gtrt'r' + 5,, (/)' + 2g,,,r'z' + .g,, {z'f , (CI) 

where A is the affine parameter along the geodesic, and a prime denotes differentiation with respect to A. Since we 
are interested in null trajectories, we set C — Q. For radial, ingoing geodesies, we have 9' ^ </)' = 0, and thus the 
geodesic equations reduce to the set: 

" P, (C2) 



A = 



Y 9rr 

2a^rg^ 
H^ 



where the dot denotes a derivative with respect to coordinate time, and H,. = dC/dr' . Then, starting at a certain 
value of r = ro, and for each grid point along the z direction, equations (|C2|) are integrated backwards in time using 
a second order Runge-Kutta scheme. 

Following null geodesies along z = const, lines does not guarantee that we are tracing the causal past of r = rg, 
though for the spacetimes studied here, and the coordinate system used, this should offer a good approximation. 
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Furthermore, although radial geodesies might not be the best choice, since the event horizon is an attractor, they 
will trace it accurately. (Fo r related discussions of approximate event horizon location in the axisymmetric four- 
dimensional case, see [29ll3?ill3lj .l 
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FIG. 4: Embedding diagrams of the apparent horizon, with the two angular dimensions 6 and suppressed, from the h/A 
evolution of a perturbed black string with L = lALc. These plots thus describe the intrinsic geometry of the apparent horizon, 
at the given instants of constant t, in a coordinate system with metric ds^ = df^ + dz^. Here, z is a periodic coordinate, and f 
is the areal radius of z = constant sections of the horizon. To better illustrate the dynamics of the horizon, we have extended 
the solution using the ^-periodicity, showing roughly two periods of the solution. See Fig. |5] for a plot of the length of one 
period of the apparent horizon versus time. 
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FIG. 5: The proper length of the apparent horizon curve in the (r, z) plane (between z = and « = L) as a function of time, 
from the h/4 evolution of a perturbed black string with L = lALc- 
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FIG. 6: The normalized Kretschmann invariant /"^ = IR\^/12 il4ti . evaluated on the apparent horizon of the perturbed 
black string spacetime with L = lALc {h/4), at the same times as shown in the embedding diagram plots (Fig. 2J- Note 
however, that here the horizontal axis is the coordinate z, and in particular the flat region of the curve between z « 3.5 and 
z ~ 6.5 in the last frame corresponds to the long, thin neck region shown in the embedding diagram plot. This demonstrates 
the rather severe "grid-stretching" problems we have then. For the static black string spacetime, J*^^ — 1 (shown for reference 
as a dotted line in the figure), while for a static 5D spherical black hole it evaluates to 6. This diagram therefore also supports 
the conclusion that at late (simulation) times the solution is tending towards a configuration describable as a sequence of black 
holes connected by thin black strings. 
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FIG. 7: Plots of the apparent horizon (labeled AH) and estimates of the event horizon location (CI, C2 and C3) in coordinate 
space (in contrast to the embedding coordinates used in Fig. |1J, from the evolution of a perturbed black string with L = lALc, 
computed at resolution h/A. Here, the CI (C2) curve marks the inward-directed past light-cone of the surface r = 10 (r = 4) 
ai t — 164. C3 denotes the outward-directed past of a surface just inside the apparent horizon at t = 164. Thus, moving 
backwards in time, these curves should asymptote towards the event horizon of the spacetime. These plots suggest that for 
most of the evolution (at least), the apparent horizon is an excellent approximation to the event horizon. 



